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A Note  on  the  Operator  Compact  Implicit  Method  for  the  Wave  Equation 


The  results  of  this  report  represent  a continuation  of  an  effort  to  develop  high 
accuracy  efficient  methods  for  solving  hyperbolic  equations.  Here  an  efficient 
fourth  order  method  for  the  multi-dimensional  wave  equation  Is  presented. 

These  techniques  are  also  being  used  to  solve  parabolic  equations  arising  in 
viscous  fluid  flow  problems  for  NAVAIR. 

This  study  has  been  supported  by  the  Naval  Surface  Weapons  Center  Independent 
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INTRODUCTION 

In  a previous  paper  [1]  a factorization  technique  which  utilized  compact 
Implicit  spatial  and  temporal  approximations  to  the  second  order  wave  equation 
was  developed.  The  method  proceeded  by  separately  Implementing  the  so  called 
compact  implicit  fourth  order  approximations  for  each  of  the  Individual  derivative 
terms.  Notwithstanding,  the  Implicit  nature  of  the  basic  approximation  Involved, 
a factorization  technique  was  described  which  allowed  one  to  resolve  higher  space 
dimension  problems  by  requiring  merely  the  solution  of  trldlagonal  equations.  In 
Section  VI  of  [1]  we  observed  a peculiar  aspect  of  our  factorization  approach. 

The  same  approach  required  twice  as  much  work  when  mixed  order  (first  and  second) 
spatial  derivative  terms  were  present. 

Upon  further  examination  of  this  approach  it  became  apparent  that  the 
problem  was  numerically  Improperly  posed  In  requiring  too  much  additional  data  to 
complete  the  factorization. 

In  this  note,  we  observe  that  by  changing  the  underlying  spatial  approximation 
when  lower  order  terms  are  present  It  la  possible  to  obtain  an  algorithm  which 
completely  resembles  our  algorithm  for  the  case  when  no  lower  order  terms  are 
present. 

In  a future  paper  [2]  the  operator  compact  Implicit  method  Is  developed.  In 
much  the  same  way  as  here,  for  parabolic  problems. 

SPATIAL  DISCRETIZATIONS 

The  classical  finite  difference  approach  for  solving  two  point  boundary 
value  problems  of  the  form 

(2.1)  L(u)  - au_  + bu  • f,  xc(0,l] 

XX  X 

with  u(0),  u(l)  given  Is  to  separately  substitute  standard  approximations  for  the 
first  and  second  derivatives  In  (2.1)  and  then  solve  the  resulting  system  of 
equations.  Accordingly,  the  fourth  order  compact  Implicit  scheme  applied  to  the 
solution  of  (2,1),  requires  that  one  solve  x 

(2.2)  Lh[0J ] - (I  + k v>"  ^ W1  + £ **r'  IT  “l  ' f J 

TT]  M.  Clment,  S.  H. Leventhal , "Higher  Order  Compact  Implicit  Schemes  for  the 
Vave  Equation,"  Math,  of  Comp.  Vol.  29,  No.  132,  1975,  pp.  985-994. 

[2]  M.  Clment,  S.  H.  Leventhal  and  B.  C.  Weinberg,  "The  Operator  Compact  Implicit 
Method  for  Parabolic  Equations,"  to  appear. 
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The  notation  here  picks  up  from  the  notation  of  [1].  References  to  an  equation  in 
[1]  will  be  superscripted  with  an  asterisk.  The  appearance  of  two  Implicit 
matrices  («ee  (2.2)*,  (6.2)*)  "trapped  inside  the  a^b^"  creates  problems  in 

trying  to  solve  the  resulting  linear  system  of  equations.  Several  people  have  tried 
block  methods  when  these  basic  tarms  appear  in  problems  [3].  A suitable  trl- 
dlegonal  relationship  however  can  be  obtained  by  merely  abandoning  attempts  to 
represent  the  separate  derivative  terms.  The  approach  we  adopt  is  to  represent 
L(Uj)  on  three  adjacent  points  up  to  highest  order  accuracy  possible  in  relation 
to  U3  on  the  same  three  points. 

A Taylor  serlas  analysis  shows  that  for  L(u)  on  a uniform  grid  fourth  order 
accuracy  can  be  obtained  by 


. r4  U4J.1  ♦ r?  “4  + r4  u4  1 

(2.3a)  q+  L(u)j+1  + qj  L(u)j  + qj  Ku)^  - 1 ^ *=* 


where 


(2.3b) 


(2.3c) 


qj  " 6*J  *J-1  + h(5*J-l  bJ  ' 2*J  bJ-l*  " h2  bj  bJ-l 

q®  - *U5aJ+1  aj_x  - b^  - bj+1  aj_1)  - h2  bJ+1  bJ_1J 

" 6*3  aj+l  " h(5aJ+l  bJ  ’ 2*J  bj+l)  " h2  bj  bJ+l 

rj  " IIqj(2aj+l  + 3h  bJ+l>  + qS(2aj  + hbj)  + qJ(2*j-l  " hbj-l)] 
rj  " I[qJ(2aj+l  + hVl>  + qj(2aj  " hbj>  + ^(2*j-l  * ^j-i^ 


The  above  relationship  can  be  expressed  in  an  operator  form  by  defining 
tridiagonal  displacement  operators  so  that  the  equation  (2.3)  is  represented 

by 


(2.4)  0 L(U).  - — R 0, 

^ j h2  * 3 

Kota  well,  for  conciseness  of  notation  we  are  using  here  to  represent  a 
different  operator  from  what  represented  in  [1].  (Indeed  in  the  case  that 

T5T  1.  Hirsh,  ^Higher  Order  Accurate  Difference  Solution  of  Fluid  Mechanics 
Problems  by  a Compact  Differencing  Technique,"  J.  Comp.  Physics,  19,  1975, 
pp.  90-109. 
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» 


a = 1,  b = 0 then  both  are  Identical.)  Below,  however,  we  will  still  retain  the 
definition  of  ■ (I  + 6^2)  from  our  previous  paper.  The  above  relationships 

(2.3)  were  first  presented  by  Swartz  [4].  The  operator  compact  Implicit  (PCI) 

representation  0 -1  R for  L(u)  presents  the  same  formal  appearance  as 
{ 2 * x 

Q"1  -2-  did  for  u alone.  By  noting  this  formal  similarity  It  la  clear  that 

'St  h2  " 

the  OCI  spatial  discretization  can  be  fully  Implemented  Into  our  factorization 
algorithm  for  the  time  dependent  second  order  wave  equation.  This  all  depends  on 
the  underlying  lnvertlbllity  of  Q^.  As  Swartz  observes,  for  the  case  that  a and 
b are  conatanta,  Is  invertible  so  long  as  the  so  called  mesh  Reynolds  number  r 

satisfies 

(2.5)  |r|  = |j*|  : m 

OCI  APPLIED  TO  THE  WAVE  EQUATION 

To  solve  (6.1)* 

(3.1)  utt  - L^(u)  + Ly(u) 

where  Initial  and  boundary  data  are  prescribed  and  where 


Lx(u)  = aux* 

L (u)  i cu 
7 77 


+ bu 

X 

+ du 


one  substltutea  the  operator  compact  implicit  approximation  (2.4)  for  the  respec- 
tive spatial  terms  and  the  same  compact  Implicit  scheme  as  before  for  the  temporal 
term  to  obtain 


(3.2)  (rt‘  ^ - 


«£)-!«» 


T n «C> 

* + -2- 

J,m  . 


n'-lR" 


X u“ 

J.l 


where 

«$>'■  Ri 
«$>■■  v ' 


[4]  B.  K.  Swartz,  "The  Construction  of  Finite  Difference  Analogs  of  Some  Finite 
Element  Schemes,"  Mathematical  Aspects  of  Finite  Elements  In  Partial 
Differential  Equations  (C.  DeBoer,  Ed.)  Academic  Preas,  1974,  pp.  279-312. 
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Following  our  approach  In  [1]  a factorization  of  (3.2)  can  be  accomplished 
by  adding  the  fourth  order  ten 


The  resulting  fon  (analogous  to  (2.9)*)  la 

0.3)  S [i  - <V>  ^ (q”*1)-1  »5Ti 

■ *=?..  - + *; + «?>•■  *”>  ”5.. 

The  numerical  Implementation  of  (3.3)  above  la  completely  analogous  to  the 
algorithm  In  [1]  for  equation  (4.3)*.  Again  with 

(3.4)  [I  - Rf1)  3 z£ 

solve  first  for  Z?**-  and  then  for  U?**. 

J »®  J »■ 

The  first  two  terms  on  the  right  hand  side  are  known  from  previous  time 
steps  and  need  not  be  computed.  The  other  two  terms  on  the  right  hand  side  are 
obtained  In  a manner  analogous  to  (4.5)*,  (4.6)*.  By  observing  that 

0.3)  x’tojr1  r;  - z”,.) 

the  left  hand  side  Is  obtained  without  much  work.  v“>n  5 12(Q^)_1  R^  u"  n Is 
obtained  by  solving  a tridiagonal  system. 


REMARKS 

1.  Hare  too,  the  algorithm  for  (3.3)  requires  that  one  generate  Initial 
data  for  Gj  , G®  . These  are  obtained  by  solving  trldlagonal  systems  as  in 
(4.9)*. 

2.  Boundary  Data:  The  computation  of  the  Intermediate  boundary  conditions 

may  ba  obtained  as  in  the  previous  paper.  However,  the  following  simplification 
Is  possible.  As  before  at  the  four  corner  points  use  the  analytic  representation 

of  Zj  ^ and  one  sided  differences  as  an  approximation  for  Zj Then  on  x - constant 
lines  solve  (3.4) 
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Observe,  Chet  the  x-sveep  part  of  the  algorithm  may  also  be  Implemented  on 

y * constant  boundary  lines  as  well  as  on  the  Interior  lines.  However,  this 

leads  to  s problem.  Experiments  reveal  that  It  Is  necessary  to  compute  z”+X  on 

J t® 

the  y ■ constant  boundaries  as  accurately  as  possible.  Thus  It  Is  necessary  to 

y|Xf 

uae  the  exact  data  for  U that  Is  given  on  these  lines.  However,  the 

n+1  J ,n 

definition  of  G,  * requires  that 
j 

x2[(QUrl  Rx  + (Qy)_1  RyJ  Uj,m» 

be  evaluated  and  A2(Qn)-J  Rn  must  be  obtained  from  (3.5),  which  includes  z" 

y y j 

which  Is  not  exact.  However,  this  problem  may  be  avoided  by  noting  that 


*2[«£rl  R"  + (Q^)'1  Ry]  - At2[Lun  + Lyun  + OCh^l 

G^1  may  then  be  redefined  on  y ■ constant  boundary  lines  as 
3 *“  


r.iH-1 


y 

12  ,wt  ~j,m 


At2[utt  + 0(h**)  ] . 


before  Implementing  the  first  step  of  the  algorithm. 

A Fourier  stability  analysis  of  (3.3)  reveals,  by  a perturbation  argument, 
that  so  long  as  r ^ h *v<  k that  again  for  C * max  (a,b)  and 

/C  A 1 7%  - 1 


the  amplification  factors  p satisfy 
Ip | S l + 0(At) 

and  so  the  scheme  Is  stable  in  the  sense  of  von  Neumann. 

Numerical  »«mple:  In  this  section  a numerical  example  Is  presented  demonstrating 

the  accuracy,  effectiveness  and  the  stability  of  the  method. 

Let  ft  be  defined  by  {0  £ x,  y £ .5] and  define  the  coefficients,  initial 
conditions,  and  boundary  conditions  of  (3.1)  by 


» 
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•(x.y.t) 

b(x,y,t) 

e(x,y,t) 

d(x,y,t) 

u(x,y,o) 


Cx»D2 

4(t+l)2 

4(t+l) 

w2 

4(t+l)2 

(y»l)2(xfl) 

4(t+l) 

#(xfl)(y+l) 


ut(x,y,o)  - (xfl)(y+l)e(x+1)(yfl) 


The  exact  eolutlon  Is 


u(x,y,t) 


(x+l)(y+l)(t+l) 

© 


The  net hod  was  run  for  a sequence  of  spetlal  meshes  and  time  steps,  so  that  In 
' each  subsequent  run  we  halved  the  mesh  size  and  time  step.  In  table  1 the 
accuracy  of  the  method  is  desionstrated  for  L^-error  and  relative  max-error. 


# Time 
Steps 

h 

k 

mi 

BUI 

Relative 

Max-Error 

Relative 

Max-Rate 

20 

.1 

.05 

3.847-05 

5.049-06 

3.74 

3.91 

40 

.05 

.025 

2.886-06 

3.341-07 

3.97 

3.98 

80 

.025 

.0125 

1.843-07 

2.118-08 

3.87 

3.86 

160 

.0125 

.00625 

1.264-08 

1.455-09 

Table  1 
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